Synergistic epistasis among cancer drivers can rescue early tumors from the accumulation of deleterious passengers

Epistasis among driver mutations is pervasive and explains relevant features of cancer, such as differential therapy response and convergence towards well-characterized molecular subtypes. Furthermore, a growing body of evidence suggests that tumor development could be hampered by the accumulation of slightly deleterious passenger mutations. In this work, we combined empirical epistasis networks, computer simulations, and mathematical models to explore how synergistic interactions among driver mutations affect cancer progression under the burden of slightly deleterious passengers. We found that epistasis plays a crucial role in tumor development by promoting the transformation of precancerous clones into rapidly growing tumors through a process that is analogous to evolutionary rescue. The triggering of epistasis-driven rescue is strongly dependent on the intensity of epistasis and could be a key rate-limiting step in many tumors, contributing to their unpredictability. As a result, central genes in cancer epistasis networks appear as key intervention targets for cancer therapy.


S1 Text: Supplementary Information
Synergistic epistasis among cancer drivers can rescue early tumors from the accumulation of deleterious passengers.

Approximate analytical study of the model
In this section we present an approximate analytical study of the model described in the main text, following the same general lines of reasoning as in the original works by McFarland et al. [1,2].The expressions derived in this section make heavy use of classical approximations from statistic physics (e.g., timescale separation) and population genetics (e.g., the fixation probability of beneficial and slightly deleterious mutations).The appropriateness and validity of these approximations is thoroughly discussed by McFarland et al. [1,2] and supported by simulations (S1 Fig) .We have attempted to provide a sufficient level of detail to allow the reader to follow the main derivations.In any case, we refer the reader to the original articles for further explanations.

General considerations
To complement the simulations described in the main text, we developed a series of approximate analytical expressions that helped us gain an intuitive understanding of the factors that determine tumor rescue in the presence of positive epistasis.To obtain such expressions, we adopted a number of approximations that facilitate analytical calculations.Most of the approximations have been previouly justified in the context of a simpler model without epistasis [1,2].However, since they concern general aspects of the modeling framework (for example, the existence of a critical size for tumor progression), they are equally applicable to a scenario with epistasis.Henceforth, we adopt the following assumptions and approximations: Despite the simulations operate in discrete generations and deal with discrete cell numbers, for our analytical calculations we assume that the dynamics of large populations can be approximated as a deterministic process with continuous time and population size.If the number of mutations remains stable, the population dynamics is described by the following differential equation where N is the population size, B is the cell duplication rate, and D(N ) is the cell death rate.
To maintain consistency with the simluations, we rescale time by the duplication rate B. This is equivalent to measuring time in units of cell generations.Because cell division rates vary across tissues and empirical mutation frequencies are better predicted by the number of divisions than by absolute time [3,4], measuring time in terms of cell generations is also a more appropriate choice from a biological perspective.The change of time scale leads to the following equation for the population dynamics: Note that the dynamics only depend on the ratio between death and duplication rates.Therefore, both rates can be jointly scaled by any constant factor without changing the dynamics.
As described in the main text, cells can harbor different numbers of passenger mutations (n p ), "regular" drivers (n d ) whose fitness effects are not affected by epistasis given the current mutational background, and "enhanced" drivers (n ϵ ) whose fitness effects are augmented as a result of synergistic interactions with other drivers.A driver mutation affecting a particular gene g is enhanced if there is at least another driver mutation affecting any of the genes connected to g in the epistasis network.We assume that mutations modify the fitness of the tumor by a multiplicative factor, such that where s d and s p are the fitness effects of regular drivers and passengers, respectively, and ϵ is a multiplicative factor that accounts for the effect of synergistic epistasis on the fitness effect of drivers.Moreover, we assume a density-dependent death rate where N 0 represents the carrying capacity in the absence of mutations.Combining eqns. 2 and 4, we obtain the following logistic equation for the population dynamics: From eqns. 3 and 5 it becomes clear that mutations affect the system by modifying its effective carrying capacity.
Based on eqn.5, the equilibrium population size (N * ) with constant mutation load is equal to Eqn. 5 also allows estimating the duration of the transient before that equilibrium size is reached.For example, starting from 1 cell, the population reaches 99.9% of the initial carrying capacity (set to N 0 = 1000) in just 14 generations.Considering the range of fitness effects that we explored in the simulations (Table 1 in the main text) and their contribution to B, when a new mutation reaches fixation in the population the new equilibrium size is typically reached in very few generations.As an upper bound, let us consider the extreme case with s d = 0.1, ϵ = 10, and suppose that the fixation of an enhanced driver in the whole population increases BN 0 from 1000 to 2000 (this is the largest possible change compatible with the range of parameters and stop criteria).According to eqn. 5, the new equilibrium (within a tolerance of 0.001%) will be reached in only 7 generations.
In the case of passengers, the transient lasts < 1 generation.In contrast, the typical waiting time for the occurrence and fixation of passengers is approximately 20 generations (with µ = 10 −8 and T p = 5 × 10 6 ) and around 10 − 100 times higher for drivers (see below for details).Because of this difference in timescales, we assume that every time that a new mutation reaches fixation, the population size instantaneously adjusts to the new equilibium value given by eqn. 6. (Note that this approximation only concerns analytical calculations: the simulations explicitly consider each cell in the population and their stochastic dynamics.) We assume that driver and passenger mutations occur during cell division.This assumtion is justified by the empirical observation that cancer risk is proportional to the number of stem cell divisions in different tissues [3].The number of drivers and passengers acquired by each daughter cell follows a Poisson distribution with mean µT d for drivers and µT p for passengers, where µ is the mutation rate per site per cell per generation, T d is the number of possible driver sites per cell, and T p is the number of possible passenger sites per cell.As indicated in Table 1 (main text), T p ≫ T d .
Driver mutations are randomly assigned to cancer genes.We assume that there are G cancer genes that equally contribute to the total number of driver sites.Only a fraction f e of the cancer genes are connected to the epistasis network.In consequence, only a fraction f e of all possible drivers can be enhanced.We study two classes of networks that contain genes that, if mutated, modify the fitness effect of drivers in any other gene of the network.Those are clique-like networks, in which every gene is connected to each other; and star-like networks, in which a single central gene (the hub) is connected to all the other (peripheral) genes.In the clique, any driver mutation will trigger the effects of epistasis, modifying the fitness effects of drivers in any other gene of the network.In contrast, in the star, only mutations in the hub can act as triggers.The existence of a single trigger event implies that we can characterize the global state of clique-and star-like epistasis networks as "active" (if the trigger has occured) or "inactive" (before the trigger occurs).Note that this is not generally possible in more complex networks.
Although the simulations consider the fully stochastic nature of the mutation process, for the analytical calculations we assume that accumulation of drivers and passengers can be approximated as deterministic processes, whose rates are given by the product of the population-level mutation rates and the fixation probabilities.This approximation is justified by numerical simulations, as shown by McFarland et al. [1,2].It follows from the previous point that the population-level mutation rates are µ(1 − f e )T d N for regular drivers, µf e T d N for enhanced drivers, and µT p N for passengers.To estimate the probability of fixation (P f ix ), we apply the classical expression for populations with constant size [5], originally derived for Wright-Fisher and Moran models.
We further simplified this expression by taking the limit of large population size (N → ∞) in the case of drivers and quasi-neutrality (s p → 0) in the case of passengers [2].
for regular drivers With these approximations, the rate of mutation accumulation becomes The different timescales of the population dynamics in the absence of new mutations (eqn.6) and the mutation accumulation process (eqn.9) allow us to approximate the dynamics of the cell population size in the presence of mutations as

Critical population size
In a previous work, McFarland et al. showed that, in the absence of epistasis, this model is characterized by the existence of a critical population size, N c .If that critical population size is not reached, the tumor recedes with very high probability [1].This phenomenon results from a tug-of-war between beneficial drivers and slightly deleterious passengers (this interpretation becomes clear from eqn. 10).If the population is sufficiently small, the rate at which new drivers appear and reach fixation is not enough to compensate for the fitness loss due to the accumulation of passengers, which leads to the decay of the population.Here we use the same approximations as in [1] to obtain the critical population size in the presence of positive epistasis and show that its value is smaller than in the absence of epistasis.To obtain the critical population size, we replace eqns.11 and 9 into eqn.10.After grouping terms, we obtain: where N ′ c can be interpreted as a critical population size (see below) and has the following expression: In the absence of epistasis (f e = 0 or ϵ = 1), the critical population size becomes As anticipated, N ′ c < N c if ϵ > 1 (positive epistasis).Notably, the solution of eqn. 12 is an increasing function of t (with lim t→∞ c and a decreasing function of t (with lim t→∞ N (t) = 0) if N (0) < N ′ c .Thus, in the deterministic version of the model described by eqn.12, the critical population size separates tumors in two classes: those that grow indefinitely (N > N ′ c ) and those that become extinct (N < N ′ c ).Although our notation does not make it explicit, it is important to note that, because in deriving eqn. 10 we assumed that the population size was always at equlibrium for a given mutational load, the variable N in eqns.10-12 corresponds to the equilibrium states of eqn. 5 (denoted by N * in eqn.6).Therefore, the dynamics of the "population size" (N ) in eqns.10-12 represents, in fact, the temporal evolution of the effective carrying capacity of the system, that takes initial value N 0 (regardless of the actual population size at the beginning of the simulation) and changes with the mutational profile as BN 0 .Consistent with this interpretation, stochasatic simulations starting from a single cell usually progress to cancer if the basal carrying capacity N 0 is greater than N ′ c (except for a small probability of extinction of order 1/N in the few first generations) [2].Similarly, simulations starting with a large number of cells usually become extinct if N 0 < N ′ c .This happens because, as discussed above, the population typically reaches its equilibrium size N * = N 0 before any drivers or passengers can reach fixation, effectively erasing any memory of the actual initial size.

Critical driver-to-passenger ratio
In this section, we derive a relationship for the ratio of drivers to passengers that separates growing from receding tumors.A similar expression was obtained by McFarland et al. [2] in the absence of epistasis.
Here, we extend it to tumors with positive epistasis.We start from the expression for the equilibrium population size, N * = BN 0 (eqn.6): Note that we have used N instead of N * for simplicity and consistency with the notation used in eqns.10-12.Now we impose the condition for tumor growth in the absence of epistasis, derived in the section 1.2: where Eqn. 16-17 define a critical line, R c , in the (n p , n d ) plane, such that, in the absence of epistasis, mutation configurations above the line lead to tumor progression and those below the line do not [2].We now intend to obtain an analogous expression that accounts for epistasis.To make it applicable to real-world scenarios in which the epistasis network may be unknown, the expression should include regular and enhanced drivers in a single variable, that we denote as n d,ϵ .To that end, we consider two extreme cases: (i) The probability of fixation of enhanced drivers is the same as regular drivers (this is true, for example, in the star before the trigger occurs).In that case, the fraction of enhanced and regular drivers is proportional to the fraction of cancer genes in and out of the epistasis network, that is n ϵ = f e n d,ϵ and n d = (1 − f e )n d,ϵ .(ii) The probability of fixation of regular drivers is negligible with respect to enhanced drivers (this would occur, for example, in a strict two-hit scenario).In that case, n ϵ = n d,ϵ and n d = 0.The critical conditions derived from these two scenarios are, respectively, In both scenarios, the factor to the left of R c is smaller than 1.Therefore, epistasis reduces the critical number of drivers relative to passengers required for cancer progression by reducing the intercept and the slope.

Tumor rescue
Driver mutations in central genes of the epistasis network (trigger drivers) can facilitate tumor rescue by changing the critical population size from N c to N ′ c as the epistasis network "switches on".To illustrate that possibility, let us consider a population with basal carrying capacity N 0 < N c .According to eqn. 12, the population will initially recede.Let us now suppose that a trigger driver reaches fixation at time t D .From that moment on, drivers hitting any gene of the epistasis network will have an enhanced fitness effect ϵs d and the dynamics of the population will be described by eqn. 12 with critical size N ′ c (given by eqn.13).Moreover, the fixation of the trigger driver produces a change ∆N D in the population that we assume to be effectively instantaneous.If the new population size fulfils the condition N (t D )+∆N D > N ′ c the tumor will expand.We term this phenomenon tumor rescue.
Unlike other drivers and passengers, trigger drivers are unique random events whose study requires a probabilistic approach.Let P (t D = t) be the probability density function for the time at which the trigger driver reaches fixation.Following the previous argument, the probability of tumor rescue can be calculated as where ) is the Heaviside step function and represents the condition for supercritical proliferation after the trigger mutation.To obtain the overall rescue probability, we need to calculate the tumor size just before the trigger driver, N (t D ), the increase in tumor size due to the fixation of the trigger driver, ∆N D and the probability that the trigger driver appears at any given time, P (t D = t).In the following subsections we will derive approximate expressions for these values.

Population size immediately before triggering
To calculate the population size at the trigger time, we reproduce the same line of reasoning that led us to derive eqns.12-14 but with a subtle modification.Because we are conditioning the dynamics to not having a trigger, the target size for drivers must be modified to exclude every mutation that could act as a network-level trigger.In the case of cliques, that implies excluding driver mutations in any gene of the network, whereas in the case of stars, that implies excluding drivers in the hub.As a result, the pre-trigger target sizes for drivers in cliques and stars become, respectively With these definitions, the dynamics of the population conditioned to the absemce of a trigger obeys the following differential equation with x ∈ {c, s} for cliques and stars, respectively, and Integrating eqn.23 from t = 0 to t = t D we obtain the expression for the population size at the trigger time: x N (t D ) = x N c 1 + ( x N c /N 0 − 1)e µTpspt D (26)

Population size immediately after triggering
The next step is to calculate the population size immediately after the fixation of the trigger driver.We denote the post-trigger population size as x N D and it is equal to the pre-trigger population size plus the increase in population size due to the trigger: with x ∈ {c, s} for cliques and stars, respectively.The expression for cliques can be readily obtained by applying eqn.11.In a clique, the trigger driver is the first driver hitting any gene from the network.In consequence, the immediate fitness effect of the trigger is just s d (it is not enhanced by epistasis because there are no mutations in neighbor genes).For that reason, we have that The case of stars is more complicated.Indeed, by the time that a trigger driver hits the hub, other drivers can already be present in peripheral cancer genes.As a result, the fitness increase associated with the appearance of the trigger driver depends on the number of peripheral drivers n ϵ .Taking that into account, the population size immediately after the trigger driver reaches fixation in stars networks becomes: where n ϵ (t D ) is the number of peripheral drivers that have reached fixation at time t D .(The first factor of the second case corresponds to the fitness effect of the trigger driver and the second factor accounts for the enhanced fitness effect of preexisting peripheral drivers.)A rigorous study of the star would require reformulating eqn.20 to consider cases with and without preexisting peripheral drivers and their respective probabilities.Instead, for simplicity, we replaced the random variable s N D (= s N (t D ) + s ∆N D ) in the Heaviside step function of eqn.20 by a weighted average ⟨ s N D ⟩ that takes into account the probabilities of having zero or more prefixated periphearal drivers.Although the error introduced by this approximation is difficult to estimate a priori, we validated this approach by comparing the resulting values of the rescue probability with those obtained thorugh simulations (S1 Fig, right).We calculated the weighted average of s N D as where P (n ϵ ) is the probability to find n ϵ prefixated drivers in peripheral cancer genes.
Because occurrence and fixation of peripheral drivers before triggering is independent of the presence of other peripheral drivers, P (n ϵ ) can be approximated by a Poisson distribution with mean λ(t D ) = µf e T d s d 1+s d t D τ =0 N (τ ).(Note that the expression of λ results from summing up the mean number of drivers per generation per cell for all cells in all generations, from t = 0 until t D .)As a result, the population size after fixation of the trigger driver, ⟨ s N D ⟩ becomes (where we have omitted the explicit dependency of λ on t D to lighten notation.)By defining A = 1+ϵs d 1+s d and operating: This expression can be simplified by taking into account that e −Aλ (Aλ) nϵ nϵ! corresponds to the probability mass function of a Poisson distribution with parameter Aλ.As a result, Finally, after some manipulation and replacing back A: Let us define the critical trigger time t c as the value of t D that solves eqn.35.Note that, if t c < 0, the tumor cannot be rescued because c N D < N ′ c for any positive trigger time.If t c > 0, we have that Therefore, the integral in eqn.20 can be rewritten as This formulation of the rescue probability highlights the fact that, in the deterministic approximation, rescue is only possible if the trigger driver occurs before t c .This observation justifies the interpretation of t c as a critical time.
Extnding these arguments to the star is not straightforward because the complex dependency between ⟨ s N D ⟩ and t D (ultimately due to the accumulation of peripheral drivers over time) does not allow rule out the possibility that the argument of the Heaviside step function becomes negative for some t D < t min .However, we have numerically verified that in all the scenarios explored in this work, the equation s N D (t D ) = N ′ c has a single solution, with the argument of the Heaviside step function changing from positive to negative as t D increases.In consequence, eqn.37 also applies in the case of the star.To apply eqn.37 to the star, t c must be first obtained by numerically solving the folowing equation for the implicit variable t D :

Rescue probability
It follows from eqn. 37 that the overall rescue probability P SOS is equal to the probability that the trigger driver reaches fixation before the critical time t c .By analogy with the arguments that led to eqn. 9, the rate at which trigger drivers appear and reach fixation is , where T D is the mutational target size associated with the trigger driver (f e T d in cliques and T d /G in stars) and s D is the fitness effect of the trigger driver (whose possible values are derived below).If we consider trigger fixation as a stochastic process with rate q, the probablity that no driver has reached fixation after a time t obeys the following differential equation dP no trigger (t) dt = −q(t)P no trigger (39) (note the dependency of q(t) with t through N (t) and, in the case of the star, s D , as described below.)By integrating this equation, we obtain a general expression for the probability that the trigger driver occurs before a given time t.
Replacing t by t c and the full expression for q we get The fitness effect s D is, by definition, the solution of the equation where B and B ′ are the birth rates before and after acquiring the trigger driver, respectively.According to eqn. 3, the fitness effect of the trigger driver in clique ( c s D ) and star ( s s D ) epistasis networks are As previously discussed in the context of eqn.29, s s D depends on the presence of drivers in peripheral genes of the star when the trigger occurs.Moreover, s s D implicitly varies with time because the probability to have n ϵ periphearal drivers, P (n ϵ ) varies with time.As we did for s N D (eqn.30), to calculate P SOS in the star we replaced the term s D 1+s D by a weighted average ⟨ s D 1+s D ⟩ that takes into account the probability to find preexistent peripheral drivers.The error introduced by this approximation is difficult to evaluate a priori, although a comparison of the final estimates of P SOS and the values obtained from simulation indicates that, within the range of parameters explored in this study, the differences are generally small (S1 Fig, right).We calculated the weighted average of

Critical rescue time for population doubling
The rescue of the tumor has been defined in terms of reaching a size greater than the post-trigger critical size N ′ c .However, it is possible that, even fulfilling this condition, the growth rate of the rescued tumor is too slow to be clinically relevant.To account for that possibility, in the simulations we only classified the outcome as tumor progression if the size doubled before t max = 30000 generations.From an analytical perspective, this condition sets a second critical time t ′ c .Rescue must occur before t ′ c so that there is enough time for the population to double its size before reaching t max .The value of t ′ c is obtained by applying eqn.26 (that describes the change in population size with time) with the population at the critical rescue time as initial condition, making it equal to 2N 0 , and numerically solving the resulting equation: Finally, the probability that the population is rescued and doubles its size before t max is calculated by applying eqns.43, 45 and 48 to the smallest of t c and t ′ c .
1+s d for regular drivers v ϵ = µf e T d N ϵs d 1+ϵs d for enhanced drivers v p = µT p + O(s p ) for passengers (9) In this expression, ∆N d , ∆N ϵ , and ∆N p represent the change in the equilibrium population size as a consequence of the fixation of a new regular driver, enhanced driver, and passenger, respectively.By combining eqns.3 and 6 and comparing the equlibrium population sizes with n d and n d + 1 regular drivers (same for enhanced drivers and passengers), we obtain that∆N d = s d Nfor regular drivers ∆N ϵ = ϵs d N for enhanced drivers ∆N p = −s p N for passengers(11)

) 1 . 4 . 3
Critical trigger time Because c N (t D ) is a decreasing function of t D (eqn.26), the post-trigger population size in the clique ( c N D in eqn.28) also decreases with t D , tending to 0 as t D → ∞.As a result, the argument of the Heaviside function in eqn.20 only changes its sign at most once, from being positive to being negative.The value of t D at which the argument of the Heaviside step function becomes negative can be obtained by solving c N D (t D ) = N ′ c (35)